二次规划是指这样的问题：找到向量 x，它能在满足特定线性约束的条件下最小化二次函数：

满足 A·x ≤ b, Aeq·x = beq, l ≤ x ≤ u。

interior-point-convex 算法执行以下步骤：

预求解/后求解

生成初始点

预测-校正

停止条件

不可行性检测

注意

算法有两种代码路径。当黑塞矩阵 H 是由双精度值组成的普通（满）矩阵时，它采用一种代码路径；当 H 是稀疏矩阵时，它采用另一种代码路径。有关稀疏数据类型的详细信息，请参阅稀疏矩阵。通常，对于非零项相对较少的大型问题，该算法在 H 指定为 sparse 时执行更快。同样，对于较小或相对稠密的问题，该算法在 H 指定为 full 时执行更快。

该算法首先尝试通过消除冗余和简化约束来简化问题。在预求解步骤中执行的任务可以包括以下内容：

检查是否有任何变量具有相等的上界和下界。如果有，请检查可行性，然后固定变量的值以去除变量。

检查是否有任何线性不等式约束只涉及一个变量。如果是，请检查可行性，然后将线性约束更改为一个边界。

检查是否有任何线性等式约束只涉及一个变量。如果是，请检查可行性，然后固定该变量的值以去除该变量。

检查是否有任何线性约束矩阵具有全零行。如果是，请检查可行性，然后删除这些行。

确定边界和线性约束是否一致。

检查是否有任何变量仅作为线性项出现在目标函数中，而不出现在任何线性约束中。如果是这样，请检查可行性和有界性，然后将变量固定为其适当的边界值。

通过添加松弛变量，将任何线性不等式约束更改为线性等式约束。

如果算法检测到不可行或无界的问题，它会停止并发出相应的退出消息。

该算法可能到达一个表示解的单个可行点。

如果算法在预求解步骤中没有检测到不可行或无界的问题，并且预求解没有产生解，则算法会继续其后续步骤。到达停止条件后，算法会重新构造原始问题，这将撤消任何预求解变换。这最后一步是后求解步骤。

有关详细信息，请参阅 Gould 和 Toint 合著的 [63]。

该算法的初始点 x0 是：

将 x0 初始化为 ones(n,1)，其中 n 是 H 中的行数。

对于同时具有上界 ub 和下界 lb 的分量，如果 x0 的分量未严格位于边界之内，则分量将设置为 (ub + lb)/2。

对于只有一个边界的分量，如有必要，请修改分量，使其严格位于边界内。

执行预测步（请参阅预测-校正），进行微小校正以实现可行性，而不是执行完整的预测-校正步。这将使初始点更靠近中心路径，且不必花费完整预测-校正步的开销。有关中心路径的详细信息，请参阅 Nocedal 和 Wright 的论著 [8]（第 397 页）。

稀疏内点凸算法和满内点凸算法的主要区别在预测-校正阶段。两种算法相似，但在某些细节上有所不同。有关基本算法说明，请参阅 Mehrotra 的 [47]。

算法首先通过将 A 和 b 乘以 -1，将线性不等式 Ax = b 形式的不等式。这不影响解，但可使其形式与一些文献中的形式保持一致。

稀疏预测-校正算法

满预测-校正算法

稀疏预测-校正算法.  与 fmincon 内点算法相似，稀疏 interior-point-convex 算法尝试找到 卡鲁什-库恩-塔克 (KKT) 条件成立的点。对于二次规划定义中所述的二次规划问题，这些条件是：

Hx+c−AeqTy−A¯Tz=0A¯x−b¯−s=0Aeqx−beq=0sizi=0, i=1,2,...,ms≥0 z≥0.

此处

A¯ 是扩展的线性不等式矩阵，其中包括以线性不等式形式编写的边界。b¯ 是对应的线性不等式向量，包括边界。

s 是将不等式约束转换为等式约束的由松弛变量组成的向量。s 的长度为 m（即线性不等式和边界的数目）。

z 是对应于 s 的拉格朗日乘数组成的向量。

y 是与等式约束相关联的拉格朗日乘数组成的向量。

该算法首先根据牛顿-拉夫逊公式计算一个预测步，然后再计算一个校正步。校正尝试以更好的方式实现非线性约束 sizi = 0。

预测步的定义：

rd，对偶残差：

rd=Hx+c−AeqTy−A¯Tz.

req，原始等式约束残差：

req=Aeqx−beq.

rineq，原始不等式约束残差，包括边界和松弛变量：

rineq=A¯x−b¯−s.

rsz，互补残差：

rsz = Sz。

S 是松弛项组成的对角矩阵，z 是拉格朗日乘数组成的列矩阵。

rc，平均互补：

rc=sTzm.

在牛顿步中，x、s、y 和 z 中的变化由下式给出：

然而，完整的牛顿步可能不可行，因为存在对 s 和 z 的正性约束。因此，quadprog 会根据需要缩短步长以保持正性。

此外，为了在内部保持“中心化”位置，该算法不尝试求解 sizi = 0，而是采用正参数 σ，并尝试求解

sizi = σrc。

quadprog 用 rsz + ΔsΔz – σrc1 取代牛顿步方程中的 rsz，其中 1 是由 1 组成的向量。此外，quadprog 对牛顿方程进行重新排序，以获得对称的、数值更稳定的方程组来计算预测步。

在计算校正的牛顿步后，算法执行更多计算，以获得更长的当前步，并准备好计算更优的后续步。这些多重校正计算可以提高性能和稳健性。关详细信息，请参阅 Gondzio 的论著 [5]。

满预测-校正算法.  满预测-校正算法不会将边界合并为线性约束，因此它有另一组对应于边界的松弛变量。该算法将下界移至零。如果某个变量只有一个边界，则算法将上界不等式反向，进而将其转换为以零为下界。

A¯ 是同时包含线性不等式和线性等式的扩展线性矩阵。b¯ 是对应的线性等式向量。A¯ 还包括用于扩展向量 x 的项，这些项使用松弛变量 s 将不等式约束转换为等式约束：

A¯x=(Aeq0AI)(x0s),

其中，x0 表示原始 x 向量。

KKT 条件是

为了找到公式 3 的解 x、松弛变量和对偶变量，该算法基本上考虑牛顿-拉夫逊步：

其中，X、V、W、T 是分别对应于向量 x、v、w 和 t 的对角矩阵。方程最右边的各个残差向量分别是：

rd，对偶残差

rp，原始残差

rub，上界残差

rvx，下界互补残差

rwt，上界互补残差

该算法以如下方式来求解 公式 4：先将其转换为以下对称矩阵形式

其中

D=H+X−1V+T−1WR=−rd−X−1rvx+T−1rwt+T−1Wrub.

D 和 R 定义中的所有矩阵的逆都很容易计算，因为矩阵是对角矩阵。

要从公式 4 推导公式 5，请注意 公式 5 的第二行与 公式 4 的第二个矩阵行相同。公式 5 的第一行来自于先求解公式 4 的最后两行以得到 Δv 和 Δw，再求解以得到 Δt。

为了求解公式 5，该算法采用 Altman 和 Gondzio·的论著 [1] 中的基本思路。该算法通过 LDL 分解求解对称方程组。正如 Vanderbei 和 Carpenter 等在论著·[2] 中指出的那样，这种分解无需进行主元消去也可保持数值稳定，因此执行速度可以很快。

在计算校正的牛顿步后，算法执行更多计算，以获得更长的当前步，并准备好计算更优的后续步。这些多重校正计算可以提高性能和稳健性。关详细信息，请参阅 Gondzio 的论著 [5]。

quadprog 满预测-校正算法与 linprog 'interior-point' 算法基本相同，但前者还包括二次项。请参阅预测-校正。

[1] Altman, Anna and J. Gondzio. Regularized symmetric indefinite systems in interior point methods for linear and quadratic optimization. Optimization Methods and Software, 1999. Available for download here.

[2] Vanderbei, R. J. and T. J. Carpenter. Symmetric indefinite systems for interior point methods. Mathematical Programming 58, 1993. pp. 1–32. Available for download here.

预测-校正算法会持续迭代，直至达到一个可行（在容差范围内满足约束）且相对步长较小的点。具体来说，定义

ρ=max(1,‖H‖,‖A¯‖,‖Aeq‖,‖c‖,‖b¯‖,‖beq‖)

当满足所有下列条件时，算法停止：

‖rp‖1+‖rub‖1≤ρTolCon‖rd‖∞≤ρTolFunrc≤TolFun,

其中

rc=maxi(min(|xivi|,|xi|,|vi|),min(|tiwi|,|ti|,|wi|)).

rc 本质上是在测量互补残差 xv 和 tw 的大小，它们均为由零组成的向量，位于同一个解处。

quadprog 在每次迭代中计算一个评价函数 φ。评价函数是可行性的测度。如果评价函数变得过大，quadprog 将停止。在本例中，quadprog 声明该问题不可行。

评价函数与问题的 KKT 条件相关 - 请参阅预测-校正。使用以下定义：

ρ=max(1,‖H‖,‖A¯‖,‖Aeq‖,‖c‖,‖b¯‖,‖beq‖)req=Aeqx−beqrineq=A¯x−b¯−srd=Hx+c+AeqTλeq+A¯Tλ¯ineqg=12xTHx+cTx−b¯Tλ¯ineq−beqTλeq.

A¯ 和 b¯ 表示扩展的线性不等式系数，其中包含表示边界的项，可用于稀疏算法。表示法 λ¯ineq 与之相似，表示线性不等式约束的拉格朗日乘数，包括边界约束。这在预测-校正中被称为 z，λeq 被称为 y。

φ 评价函数是

1ρ(max(‖req‖∞,‖rineq‖∞,‖rd‖∞)+g).

如果此评价函数变得太大，则 quadprog 会声明该问题不可行并停止，同时返回退出标志 -2。

Optimization Toolbox™ 求解器中使用的许多方法都基于信赖域，这是一个简单而功能强大的优化概念。

要理解信赖域优化方法，请考虑无约束最小化问题，最小化 f(x)，该函数接受向量参量并返回标量。假设您现在位于 n 维空间中的点 x 处，并且您要寻求改进，即移至函数值较低的点。基本思路是用较简单的函数 q 来逼近 f，该函数需能充分反映函数 f 在点 x 的邻域 N 中的行为。此邻域是信赖域。试探步 s 是通过在 N 上进行最小化（或近似最小化）来计算的。以下是信赖域子问题，

如果 f(x + s) < f(x)，当前点更新为 x + s；否则，当前点保持不变，信赖域 N 缩小，算法再次计算试探步。

在定义特定信赖域方法以最小化 f(x) 的过程中，关键问题是如何选择和计算逼近 q（在当前点 x 上定义）、如何选择和修改信赖域 N，以及如何准确求解信赖域子问题。本节重点讨论无约束问题。后面的章节讨论由于变量约束的存在而带来的额外复杂性。

在标准信赖域方法 ([48]) 中，二次逼近 q 由 F 在 x 处的泰勒逼近的前两项定义；邻域 N 通常是球形或椭圆形。以数学语言表述，信赖域子问题通常写作

其中，g 是 f 在当前点 x 处的梯度，H 是黑塞矩阵（二阶导数的对称矩阵），D 是对角缩放矩阵，Δ 是正标量，‖ . ‖ 是 2-范数。存在求解公式 7 的好算法（请参阅[48]）；此类算法通常涉及计算 H 的所有特征值，并将牛顿法应用于以下久期方程

1Δ−1‖s‖=0.

此类算法提供公式 7 的精确解。但是，它们要耗费与 H 的几个分解成比例的时间。因此，对于大型问题，需要一种不同方法。文献（[42] 和 [50]）中提出了几种基于公式 7 的逼近和启发式方法建议。Optimization Toolbox 求解器采用的逼近方法是将信赖域子问题限制在二维子空间 S 内（[39] 和 [42]）。一旦计算出子空间 S，即使需要完整的特征值/特征向量信息，求解公式 7 的工作量也不大（因为在子空间中，问题只是二维的）。现在的主要工作已转移到子空间的确定上。

二维子空间 S 是借助下述预条件共轭梯度法确定的。求解器将 S 定义为由 s1 和 s2 确定的线性空间，其中 s1 是梯度 g 的方向，s2 是近似牛顿方向，即下式的解

或是负曲率的方向，

如果 gi < 0 且 ui = ∞，则 vi = –1

如果 gi ≥ 0 且 li = –∞，则 vi = 1

非线性系统公式 13 不是处处可微的。不可微分性发生在 vi = 0 时。您可以通过保持严格的可行性（例如限制 l < x < u）来避免这些不可微分的点。

公式 13 给出的非线性方程组的缩放的修正牛顿步 sk 定义为以下线性方程组的解

在第 k 次迭代中的解，其中

并且

在此处，Jv 表示 |v| 的雅可比。对角矩阵 Jv 的每个对角分量等于 0、–1 或 1。如果 l 和 u 的所有分量均为有限的，则 Jv = diag(sign(g))。在满足 gi = 0 的点上，vi 可能不可微分。Jiiv=0 在此类点上定义。这种类型的不可微分性不值得关注，因为对于这种分量，vi 采用哪个值并不重要。此外，|vi| 在此点上仍是不连续的，但函数 |vi|·gi 是连续的。

其次，反射用于增大步长。单个反射步定义如下。给定与边界约束相交的步 p，假设与 p 相交的第一个边界约束是第 i 个边界约束（可以是第 i 个上界，也可以是第 i 个下界），则反射步 pR = p，但在第 i 个分量中除外，其中 pRi = –pi。

完成预求解步骤后，active-set 算法分两个阶段进行。

阶段 1 - 获得一个在所有约束下可行的点。

阶段 2 - 以迭代方式降低目标函数，每次迭代均保持一组活动约束并保持可行性。

active-set 策略（也称为投影方法）类似于吉尔等人在 [18] 和 [17] 中所述的策略。

该算法首先尝试通过消除冗余和简化约束来简化问题。在预求解步骤中执行的任务可以包括以下内容：

检查是否有任何变量具有相等的上界和下界。如果有，请检查可行性，然后固定变量的值以去除变量。

检查是否有任何线性不等式约束只涉及一个变量。如果是，请检查可行性，然后将线性约束更改为一个边界。

检查是否有任何线性等式约束只涉及一个变量。如果是，请检查可行性，然后固定该变量的值以去除该变量。

检查是否有任何线性约束矩阵具有全零行。如果是，请检查可行性，然后删除这些行。

确定边界和线性约束是否一致。

检查是否有任何变量仅作为线性项出现在目标函数中，而不出现在任何线性约束中。如果是这样，请检查可行性和有界性，然后将变量固定为其适当的边界值。

通过添加松弛变量，将任何线性不等式约束更改为线性等式约束。

如果算法检测到不可行或无界的问题，它会停止并发出相应的退出消息。

该算法可能到达一个表示解的单个可行点。

如果算法在预求解步骤中没有检测到不可行或无界的问题，并且预求解没有产生解，则算法会继续其后续步骤。到达停止条件后，算法会重新构造原始问题，这将撤消任何预求解变换。这最后一步是后求解步骤。

有关详细信息，请参阅 Gould 和 Toint 合著的 [63]。

在阶段 1 中，算法尝试找到一个点 x，它满足所有约束而不考虑目标函数。仅当提供的初始点 x0 不可行时，quadprog 才运行阶段 1 算法。

首先，算法尝试找到一个在所有等式约束下都可行的点，例如 X = Aeq\beq。如果方程 Aeq*x = beq 没有解 x，则算法停止。如果有解 X，下一步是满足边界和线性不等式。在没有等式约束的情况下，设置 X = x0，即初始点。

从 X 开始，算法计算 M = max(A*X – b, X - ub, lb – X)。如果 M options.ConstraintTolerance，算法为辅助线性规划问题引入非负松弛变量 γ

minx,γγ

使得

Ax−γ≤bAeq x=beqlb−γ≤x≤ub+γγ≥−ρ.

此处，ρ 是 ConstraintTolerance 选项乘以 A 和 Aeq 中最大元素的绝对值。如果算法找到 γ = 0，并且一个点 x 满足方程和不等式，则 x 是一个可行的阶段 1 点。如果在 γ = 0 的情况下，辅助线性规划问题没有解 x，则阶段 1 问题不可行。

为了求解辅助线性规划问题，算法设置 γ0 = M + 1，设置 x0 = X，并将活动集初始化为固定变量（如果有）和所有等式约束。算法将线性规划变量 p 重新表示为 x 相对于当前点 x0 的偏移量，即 x = x0 + p。该算法使用阶段 2 中求解二次规划问题相同的迭代来求解线性规划问题，并使用适当的修正黑塞矩阵。

使用变量 d，问题可表示为

此处，Ai 指 m×n 矩阵 A 的第 i 行。

在阶段 2 中，活动集 A¯k，它是在解点处的活动约束（约束边界上的约束）的估计值。

算法在每一迭代 k 中更新 A¯k，以构成搜索方向 dk 的基础。等式约束始终维持在活动集 A¯k 中。计算搜索方向 dk 以最小化目标函数，同时确保该方向位于活动约束边界上。算法基于 Zk 形成 dk 的可行子空间，Zk 的列与活动集 A¯k 的估计值正交（即 A¯kZk=0）。这就确保了搜索方向（由 Zk 的列的任意组合经线性求和形成）始终位于活动约束的边界上。

算法基于矩阵 A¯kT 的 QR 分解的最后 n – l 个列形成矩阵 Zk，其中 l 是活动约束的数目且 l < n。也就是说，Zk 由下式给出

其中

QTA¯kT=[R0].

在找到 Zk 后，算法寻找新搜索方向 dk 以最小化 q(d)，其中 dk 在活动约束的零空间中。也就是说，dk 是 Zk 的列的线性组合：对于某个向量 p，d^k=Zkp。

通过代换 dk，将二次函数视为 p 的函数，可得

求此方程关于 p 的微分，可得

∇q(p) 称为二次函数的投影梯度，因为它是投影到 Zk 定义的子空间中的梯度。项 ZkTHZk 称为投影黑塞矩阵。假设投影黑塞矩阵 ZkTHZk 是半正定矩阵，在由 Zk 定义的子空间中函数 q(p) 的最小值出现在 ∇q(p) = 0 时，它是以下线性方程组的解

然后，算法执行以下求解步

xk+1=xk+αdk,

其中

dk=Zkp.

由于目标函数的二次性质，在每次迭代中，步长 α 只有两种选择。沿 dk 的单位步是 A¯k 的零空间内到函数最小值的精确步。如果算法可以采取这样的步而不违反约束，则此步就是二次规划的解（公式 18）。否则，沿 dk 到最近约束的步小于单位步，并且下一次迭代时，算法将在活动集中加入一个新的约束。沿任一方向 dk 到约束边界的距离由下式给出

α=mini∈{1,...,m}{−(Aixk−bi)Aidk},

该定义针对不在活动集中的约束，其中的方向 dk 朝向约束边界，即 Aidk>0, i=1,...,m。

如果活动集包含 n 个独立约束，且最小值的位置未知，则算法会计算满足以下非奇异线性方程组的拉格朗日乘数 λk

如果 λk 的所有元素均为非负，则 xk 是二次规划问题 公式 1 的最优解。但是，如果 λk 的任一分量为负，并且该分量不对应于等式约束，则最小化未完成。算法会删除对应于最小的负乘数的元素，并开始一次新迭代。

有时，即使求解器用尽了所有非负拉格朗日乘数，一阶最优性测度仍高于容差，或仍不满足约束容差。在这些情况下，求解器会尝试按照 [1] 中所述的重启过程来获得更好的解。在此过程中，求解器会放弃当前活动约束集，然后稍微放松约束，并构造新的活动约束集，同时尝试以避免循环（重复返回到相同状态）的方式求解问题。如有必要，求解器可以多次执行重启过程。

注意

不要尝试通过为 ConstraintTolerance 和 OptimalityTolerance 选项设置较大的值来提前停止算法。通常，求解器在进行迭代时不检查这些值，直到到达一个潜在的停止点时才检查是否满足容差。

有时，active-set 算法可能难以检测问题何时是无界的。如果无界性 v 的方向是二次项 v'Hv = 0 的方向，则可能出现此问题。为了数值稳定性和实现乔列斯基分解，active-set 算法对二次目标加上一个值很小的严格凸项。这个值很小的项使得目标函数的边界不会为 –Inf。在这种情况下，active-set 算法会达到一个迭代限制，而不是报告解是无界的。换句话说，该算法停止时会返回退出标志 0（而不是 –3）。

[1] Gill, P. E., W. Murray, M. A. Saunders, and M. H. Wright. A practical anti-cycling procedure for linearly constrained optimization. Math. Programming 45 (1), August 1989, pp. 437–474.

当您以热启动对象作为起点运行 quadprog 或 lsqlin 'active-set' 算法时，求解器会尝试跳过第 1 阶段和第 2 阶段的许多步骤。热启动对象包含处于活动状态的约束集，此活动约束集对于新问题可能是正确的或接近正确。因此，求解器可以避免运行一些迭代而直接向活动约束集添加约束。此外，初始点可能接近新问题的解。有关详细信息，请参阅 optimwarmstart。